function sf = solve_forward(param,fp,fp_parfor,h)

alpha = param(1:9);
TFP = param([10:11 end]);
gamma = param(15:19);

sim_data = NaN( 3*fp_parfor.n_simulation , 13 ) ;
sim_data_network = NaN( 3*fp_parfor.n_simulation , 5 ) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulate the model;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

grid =  cell(fp.periods);
grid_peers =  cell(fp.periods);
grid_inv0   = cell(fp.periods);
grid_inv1   = cell(fp.periods);
hc = cell((fp.periods-1),1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
% Draw Initial Conditions;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

init_draw_log_skills = fp_parfor.init_draws;
hc{1}= exp(init_draw_log_skills);

Peers = NaN(size(hc{1},1),fp.periods);
distrib_skills   = NaN(size(hc{1},1),fp.periods);
distrib_PS = NaN(size(hc{1},1),fp.periods-1);
Skills = NaN(size(hc{1},1),fp.periods);
Peers_log = NaN(size(hc{1},1),fp.periods);

Skills(:,1) = hc{1};

%Parental Education (College \in {0,1});
College_parent = fp_parfor.College;


for t=1:(fp.periods-1)

    %Random graph generation;

    [H1, H2] = ndgrid(hc{t}, hc{t});

    if t==1

        U_links1 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2);
        U_links2 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) ;

    else

        %t-1 PS stored in the loop;     
        [PS1, PS2] = ndgrid(PS, PS);


        U_links1 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H1) + gamma(3).*log(H2) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
            gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H1) - log(H2)>0).*PS1;
        U_links2 = gamma(1).*ones(size(H1)) +  ...
            gamma(2).*log(H2) + gamma(3).*log(H1) + gamma(4).*(log(H1) - log(H2)).^(2) + ...
            gamma(5).*( (log(H1) - log(H2)).^(2) ).*(log(H2) - log(H1)>=0).*PS2;

    end


    Prob_Links = (1./(1+1./exp(U_links1))).*(1./(1+1./exp(U_links2))) ;


    ADJ = ( Prob_Links >= fp_parfor.peers_forward(:,:,1,t));
    ADJ(logical(eye(size(ADJ)))) = 0;

    peers_skills = hc{t} ;

    peers =  repmat(peers_skills',[size(Prob_Links,1) ,1]) ;
    temp = ADJ.*peers;
    H_peers =  nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;
    Peers(:,t) = H_peers;

    %Compute the mean-log peers for the indirect inference moments;
    temp = ADJ.*log(peers);
    H_peers_logs =  nansum(temp,2)./sum(temp~=0 &  isnan(temp)==0,2) ;
    Peers_log(:,t) = H_peers_logs;

    %Number of Friends
    n_friends = sum( ADJ , 2 );


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
    % Simulate Parental Optimal Behavior;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

    State_vars = [hc{t} ,  H_peers ];
    value0_non_college = Interp_function_simulation(fp_parfor.estim_value0{t,1},State_vars);
    value1_non_college = Interp_function_simulation(fp_parfor.estim_value1{t,1},State_vars);

    value0_college = Interp_function_simulation(fp_parfor.estim_value0{t,2},State_vars);
    value1_college = Interp_function_simulation(fp_parfor.estim_value1{t,2},State_vars);

    %Probability of authoritarian parenting by Parental Education (non-colleve vs. college);
    PS_prob_non_college = (1./(1+1./exp(value1_non_college-value0_non_college)));
    PS_prob_college = (1./(1+1./exp(value1_college-value0_college)));

    PS_non_college = NaN(size(value1_non_college));
    PS_college = NaN(size(value1_college));

    %Simulate the behavior given random numbers ( fp_parfor.PS_forward );
    PS_non_college(isnan(PS_prob_non_college)==0) = ( PS_prob_non_college(isnan(PS_prob_non_college)==0)>= fp_parfor.PS_forward(isnan(PS_prob_non_college)==0,1,t) ) ;
    PS_college(isnan(PS_prob_college)==0) = ( PS_prob_college(isnan(PS_prob_college)==0)>= fp_parfor.PS_forward(isnan(PS_prob_college)==0,1,t) ) ;

    %Aggregate distribution of parenting style (aggregating over types);
    PS = (College_parent==0).*PS_non_college + (College_parent==1).*PS_college;

    %P=0 & No-College
    inv0{1}  = Interp_function_simulation(fp_parfor.estim_policy0{t,1},State_vars);
    inv0{1}(inv0{1}<fp.Min_Time)=fp.Min_Time;
    inv0{1}(inv0{1}>fp.Max_Time)=fp.Max_Time;

    %P=0 & College
    inv0{2}  = Interp_function_simulation(fp_parfor.estim_policy0{t,2},State_vars);
    inv0{2}(inv0{2}<fp.Min_Time)=fp.Min_Time;
    inv0{2}(inv0{2}>fp.Max_Time)=fp.Max_Time;


    %P=1 & No-College Parent
    inv1{1}  = Interp_function_simulation(fp_parfor.estim_policy1{t,1},State_vars);
    inv1{1}(inv1{1}<fp.Min_Time)=fp.Min_Time;
    inv1{1}(inv1{1}>fp.Max_Time)=fp.Max_Time;

    %P=1 & College
    inv1{2}  = Interp_function_simulation(fp_parfor.estim_policy1{t,2},State_vars);
    inv1{2}(inv1{2}<fp.Min_Time)=fp.Min_Time;
    inv1{2}(inv1{2}>fp.Max_Time)=fp.Max_Time;

    %Optimal parental investments;
    Inv = (College_parent==0).*(inv0{1}.*(PS_non_college==0) + inv1{1}.*(PS_non_college==1)) + (College_parent==1).*(inv0{2}.*(PS_college==0) + inv1{2}.*(PS_college==1));
    %Evolution of a child's skills given optimal parental behavior;
    Skills(:,t+1) = technology(t,alpha, hc{t}, Inv, PS, H_peers, College_parent,  TFP );
    %Store distribution of skills;
    hc{t+1} =  Skills(:,t+1);

    %Updating the grid points of the model when solving backward;
    grid{t} = prctile( hc{t} , fp.percentiles_kid ) ;
    grid_peers{t} = prctile( Peers(:,t) , fp.percentiles_kid ) ;
    grid_inv0{t}   = prctile( [inv0{1} ; inv0{2} ] , fp.percentiles_inv ) ;
    grid_inv1{t}  = prctile( [inv1{1} ; inv1{2} ] , fp.percentiles_inv ) ;
    
    %Store the aggregate state variables (distribution of skills and %parenting style in the economy);
    distrib_skills(:,t+1) =  Skills(:,t+1) ;
    distrib_PS(:,t) = PS;


    %%%%%%%%%%%%%%%%%%%%%%%%%;
    %Store Simulated Data;
    %%%%%%%%%%%%%%%%%%%%%%%%%;    

    ind_obs = 1+(t-1)*fp_parfor.n_simulation:t*fp_parfor.n_simulation;

    sim_data( ind_obs , 1 ) = 1:fp_parfor.n_simulation;
    sim_data( ind_obs , 2 ) = repmat(fp.ages(t),length(ind_obs),1);
    sim_data( ind_obs , 3 ) = h;
    sim_data( ind_obs , 4 ) = hc{t};
    sim_data( ind_obs , 5 ) = hc{t+1};
    sim_data( ind_obs , 6 ) = Inv;
    sim_data( ind_obs , 7)  = Peers(:,t) ;
    sim_data( ind_obs , 8 ) = Peers_log(:,t);
    sim_data( ind_obs , 9 ) = PS;
    sim_data( ind_obs , 13) = College_parent;



    sim_data_network( ind_obs , 1 ) = 1:fp_parfor.n_simulation;
    sim_data_network( ind_obs , 2 ) = repmat(fp.ages(t),length(ind_obs),1);
    sim_data_network( ind_obs , 3 ) = h;
    sim_data_network( ind_obs , 4 ) = hc{t} ;
    sim_data_network( ind_obs , 5 ) = n_friends;

end %t

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Add next-period peers into the simulated dataset;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
for t=1:(fp.periods-2)

    ind_obs = 1+(t-1)*fp_parfor.n_simulation:t*fp_parfor.n_simulation;
    sim_data( ind_obs , 10 ) = Peers(:,t+1);
    sim_data( ind_obs , 11 ) = Peers_log(:,t+1);

end

sim_data( : , 12 ) =  fp_parfor.j ;


sf.distrib_skills = distrib_skills;
sf.distrib_PS = distrib_PS;
sf.grid = grid;
sf.grid_peers = grid_peers;
sf.grid_inv0 = grid_inv0;
sf.grid_inv1 = grid_inv1;
sf.sim_data = sim_data;
sf.sim_data_network = sim_data_network;

